Fitting models to surface#
In his series of papers studying and comparing various magnetopause empirical models, [Nguyen et al., 2022] defines numerically these models in a publicly available script. We will use some of these functions to understand whether they could describe the simulation’s magnetopause. We will be initializing each of the following fits with the simulation parameters as indicated in Initial OpenGGCM simulation parameters.. For the parameter estimation we will use directly the built-in curve-fit function of SciPy.
Here we try fitting the Shue model to the full 3D data of the maximum emissivity surface:
Show code cell source
theta_grid, phi_grid = Theta[:,0,:], Phi[:,0,:]
gamma = 0
P_dyn = 12.5
P_m = 1.5
B_z = 5
Omega = np.radians(0)
phi_flat = np.ravel(phi_grid)
theta_flat = np.ravel(theta_grid)
r_flat = np.ravel(r_surface)
# Mask to keep only finite values
mask = np.isfinite(theta_flat) & np.isfinite(r_flat)
theta_flat_valid = theta_flat[mask]
r_flat_valid = r_flat[mask]
popt, pcov = curve_fit(model.MP_Shue1998, theta_flat_valid, r_flat_valid, p0=[12.5, 5])
perr = np.sqrt(np.diag(pcov))
param_names = ['Pd', 'Bz']
for name, value, error in zip(param_names, popt, perr):
print(f"{name} = {value:.4f} ± {error:.4f}")
DP = popt[0]
Bz = popt[1]
r0 = (10.22+1.29*np.tanh(0.184*(Bz+8.14)))*DP**(-1./6.6)
a = (0.58-0.007*Bz)*(1+0.024*np.log(DP))
print(f"r0={r0}, a={a}")
from sklearn.metrics import mean_squared_error
# Compute Shue model surface with fitted parameters
r_shue_fit = model.MP_Shue1998(theta_grid, popt[0], popt[1])
# Compute RMSE (mask out invalid values if needed)
mask = ~np.isnan(r_surface) & ~np.isnan(r_shue_fit)
rmse = np.sqrt(mean_squared_error(r_surface[mask], r_shue_fit[mask]))
print(f"RMSE between r_surface and fitted Shue surface: {rmse} RE")
Pd = 11.3580 ± 0.0127
Bz = 6.5704 ± 0.1762
r0=7.956991679798463, a=0.5651498286559324
RMSE between r_surface and fitted Shue surface: 0.5972698393048848 RE
Show code cell source
theta = np.arange(0, np.pi/2, 0.01)
phi = np.arange(0, 2*np.pi, 0.01)
P_dyn = popt[0]
B_z = popt[1]
r_shue = model.MP_Shue1998(theta_grid, P_dyn, B_z)
x, y, z = model.get_cartesian(r_shue, theta_grid, phi_grid)
fig = go.Figure(data=[go.Surface(
x=x,
y=y,
z=z,
surfacecolor=r_shue,
colorscale='thermal',
cmin=0,
cmax=np.nanmax(r_shue),
colorbar=dict(title="r value")
)])
fig.update_layout(
title="Shue model",
template="plotly_dark",
scene=dict(
xaxis_title="X",
yaxis_title="Y",
zaxis_title="Z"
)
)
fig.show()